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Abstract 

We develop a completely new and straightforward method for simulating the joint 
law of the position and running maximum at a fixed time of a general Levy process 
with a view to application in insurance and financial mathematics. Although differ- 
ent, our method takes lessons from Carr's so-called 'Canadization' technique as well 
as Doney's method of stochastic bounds for Levy processes; see Carr [7] and Doney 
[10]. We rely fundamentally on the Wiener-Hopf decomposition for Levy processes 
as well as taking advantage of recent developments in factorisation techniques of the 
latter theory due to Vigon [23] and Kuznetsov [14]. We illustrate our Wiener-Hopf 
Monte-Carlo method on a number of different processes, including a new family of 
Levy processes called hypergeometric Levy processes. Moreover, we illustrate the 
robustness of working with a Wiener-Hopf decomposition with two extensions. The 
first extension shows that if one can successfully simulate for a given Levy processes 
then one can successfully simulate for any independent sum of the latter process 
and a compound Poisson process. The second extension illustrates how one may 
produce a straightforward approximation for simulating the two sided exit problem. 
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1 Introduction 



Let us suppose that X = {X t : t > 0} is a general Levy process with law P and Levy 
measure II. That is to say, X is a Markov process with paths that are right continu- 
ous with left limits such that the increments are stationary and independent and whose 
characteristic function at each time t is given by the Levy-Khinchine representation 

E[e i6Xt ] = e- me) , 9ER, (1) 

where 

(0) = i6a + l -a 2 6 2 + J (1 - e Wx + i9xl {]x]<1} )U(dx) . (2) 

We have a G K, a 2 > and II is a measure supported on K with II ({0}) = and 
J R (x 2 A l)n(dx) < oo. Starting with the early work of Madan and Seneta [19], Levy 
processes have played a central role in the theory of financial mathematics and statistics 
(see for example the books [4, 8, 20, 21]). More recently they have been extensively used in 
modern insurance risk theory (see for example Kluppelberg et al. [13], Song and Vondracek 
[22]). The basic idea in financial mathematics and statistics is that the log of a stock price 
or risky asset follows the dynamics of a Levy process whilst in insurance mathematics, it 
is the Levy process itself which models the surplus wealth of an insurance company until 
ruin. There are also extensive applications of Levy processes in queuing theory, genetics 
and mathematical biology as well as through their appearance in the theory of stochastic 
differential equations. 

In both financial and insurance settings, a key quantity of generic interest is the joint 
law of the current position and the running maximum of a Levy process at a fixed time 
if not the individual marginals associated with the latter bivarite law. For example, if 
we define X t = sup s<t X s then the pricing of barrier options boil down to evaluating 
expectations of the form E[/(x + X t )l^ x+ x t>b ^\ for some appropriate function f(x) and 
threshold b > 0. Indeed if f(x) = (K — e x ) + then the latter expectation is related to the 
value of an 'up-and-in' put. In credit risk one is predominantly interested in the quantity 
P(X t < x) as a function in x and t, where P is the law of the dual process —X. Indeed 
it is as a functional of the latter probabilities that the price of a credit default swap 
is computed; see for example the recent book of Schoutens and Cariboni [21]. One is 
similarly interested in P(X t > x) in ruin theory as these probabilities are also equivalent 
to the finite-time ruin probabilities. 

One obvious way to do Monte-Carlo simulation of expectations involving the joint law 
of (X t ,X t ) that takes advantage of the stationary and independent increments of Levy 
processes is to take a random walk approximation to the Levy process, simulate multiple 
paths, taking care to record the maximum for each run. Whilst one is able to set things 
up in this way so that one samples exactly from the distribution of X t , the law of the 
maximum of the underlying random walk will not agree with the law of X t . 
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Taking account of the fact that all Levy processes respect a fundamental path de- 
composition known as the Wiener-Hopf factorisation, it turns out there is another very 
straightforward way to perform Monte-Carlo simulations for expectations involving the 
joint law of (X t ,X t ) which we introduce in this paper. Our method allows for exact 
sampling from the law of (X g ,X g ) where g is a random time whose distribution can be 
concentrated arbitrarily close around t. 

There are several advantages of the technique we use which are discussed in detail in 
the subsequent sections of this paper. Firstly, when it is taken in context with very recent 
developments in Wiener-Hopf theory for Levy processes, for example recent advances in 
the theory of scale functions for spectrally negative processes (see Kyprianou et al. [17]), 
new complex analytical techniques due to Kuznetsov [14] and Vigon's theory of philan- 
thropy (see [23]), one may quickly progress the algorithm to quite straightforward numer- 
ical work. Secondly, our Wiener-Hopf method takes advantage of a similar feature found 
in the, now classical, 'Canadization' method of Carr [7] for numerical evaluation of opti- 
mal stopping problems. The latter is generally acknowledged as being more efficient than 
appealing to classical random walk approximation Monte-Carlo methods. Indeed, later in 
this paper, we present our numerical findings with some indication of performance against 
the method of random walk approximation. In this case, our Wiener-Hopf method appears 
to be extremely effective. Thirdly, in principle, our method handles better the phenomena 
of discontinuities which can occur with functionals of the form E[/(aj + X t )l{ x+ x t>b y] at 
the boundary point x = b. It is now well understood that the issue of regularity of the 
upper and lower half line for the underlying Levy process (see Chapter 6 of [15] for a 
definition) is responsible the appearance of a discontinuity at x = b in such functions 
(cf. [1]). The nature of our Wiener-Hopf method naturally builds the distributional atom 
which is responsible for this discontinuity into the simulations. 

Additional advantages to the method we propose include its simplicity with regard to 
numerical implementation. Moreover, as we shall also see in Section 4 of this paper, the 
natural probabilistic structure pertaining to Levy processes that lies behind our so-called 
Wiener-Hopf Monte-Carlo method also allows for additional creativity when addressing 
some of the deficiencies of the method itself. 

2 Wiener-Hopf Monte-Carlo simulation technique 

The basis of the algorithm is the following simple observation which was pioneered by Carr 
[7] and subsequently used in several contexts within mathematical finance for producing 
approximate solutions to free boundary value problems that appear as a result of optimal 
stopping problems that characterise the value of an American-type option. 

Suppose that ei, e 2 , • • • are a sequence of i.i.d. exponentially distributed random vari- 
ables with unit mean. Suppose they are all defined on a common product space with 
product law P which is orthogonal to the probability space on which the Levy process X 
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is defined. For all t > 0, we know from the Strong Law of Large Numbers that 

n t 

— ej — > t as n t oo (3) 
i=i n 

P-almost surely. The random variable on the left hand side above is equal in law to a 
Gamma random variable with parameters n and n/t. Henceforth we write it g(n, n/t). Re- 
call that (X, P) is our notation for a general Levy process. Then writing X t = sup s<t X s we 
argue the case that, for sufficiently large n, a suitable approximation to ¥(X t G da;, X t e 
dy) is P x F(X g{rhn/t) e dx, X g{n ^ /t) G dy). 

This approximation gains practical value in the context of Monte-Carlo simulation 
when we take advantage of the fundamental path decomposition that applies to all Levy 
processes over exponential time periods known as the Wiener-Hopf factorisation. 

Theorem 1. For all n > 1 and A > define g(n, A) := Y^=i e «M- Then 

(X g(n , A) ,X g(niA) ) £ (V(n,X),J(n,X)) (4) 

where V(n, A) and J(n, A) are defined iteratively for n > 1 as 

V(n,X) = V(n-l,X)+S^+I^ 
J(n,X) = max(j(n-l,A) ,V(n - 1, X) + 



and V(0, A) = J(0, A) = 0. Here, = 1^ = 0, {S^ : j > 1} are an i.i.d. sequence 
of random variables with common distribution equal to that of X ei /\ and {/j^ : j > 1} 
are another i.i.d. sequence of random variables with common distribution equal to that of 

=Xei/A- 

Proof. Suppose we define X s ^ t = sup s<u<t X u . Then it is trivial to note that 

n 

X g(n,X) = \/ -^g(i-l,A),g(i,A) (5) 
i=l 

where g(0, A) := 0. 

Next we prove by induction that for each k — 0, 1, • • • 

(X sM :n<k)± {j^Sf + I?} :n<k^. (6) 

Note first that the above equality is trivially true when k = 1 on account of the Wiener- 
Hopf factorisation. Indeed the latter tells us that X ei /\ and X ei /\ — X ei /\ are independent 
and the second of the pair is equal in distribution to X ei , x . Now suppose that (6) is true for 
k > 1. Then stationary and independent increments of X together with the Wiener-Hopf 
factorisation imply that 

Ag(fc+1,A) - *g(k,\) + A e fc+1 /A - A g(fc,A) + ^A + J A 
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where X^ fc+1 ^ is an independent copy of X, S^ +1 ' := sup s<efe+i / A xi k+1 ^ and J^ fc+1 - ) 
inf s < efc+1 /A X^ +1 . The induction hypothesis thus holds for k + 1. 

For n = 0, 1, • • • , stationary and independent increments of X allows us to write 

-Xg(n,A),g(n+l,A) = -Xg(n,A) + SU P = Xg{n,X) + S^ +1 . 

s<e n+ i/X 

Hence for k = 0,1, ■ ■ ■ 

(n 
Y^iS? + I®} + S«*V : n < k 
3=1 

^From (5) and (6) the result now follows. □ 

Note that the idea of embedding a random walk into the path of a Levy process with 
two types of step distribution determined by the Wiener- Hopf factorisation has been used 
in a different, and more theoretical context, by Doney [10]. 

Given (3) it is clear that the pair (V(n,n/t), J{n,n/t)) converges in distribution to 
(X t ,X t ). This suggests that we need only to be able to simulate i.i.d. copies of the dis- 
tributions of S n /t := S^l and I n / t := I^) t and then by a simple functional transformation 
we may produce a realisation of random variables (X g (n, n /t), -Xg(n,n/t))- Given a suitably 
nice function F, using standard Monte-Carlo methods one estimates for large k 

k 

E[F(X t ,X t )\ ~ ^F(^(n,n/t), J^\n,n/t)) (7) 

m=l 

where (V^ m \n, n/t), J^ m \n, n/t)) are i.i.d. copies of (V(n,n/t), J(n,n/t)). Indeed the 
strong law of large numbers implies that the right hand side above converges almost 
surely as k t oo to E x E(F(X s ^ n)n M,X s ^ n)n / t ))) which in turn converges as n | oo to 
E(F(X t ,X t )). 

The Central Limit Theorem indicates that the right hand side of (7) converges to 
E[F(V(n,n/t), J{n,n/t))] at rate 0(k~»). Below we give an indiction of the rate of con- 
vergence of E[F{V(n, n/t), J{n,n/t))] to the desired expectation E[F(X t , X t )]. 

Lemma 1. Assume that function F(x,y) : 1R 2 i— > K belongs to the domain of C, the 
infinitesimal generator of (X t , X t ). Then as n | oo 

E[F(V(n,n/t),J(n,n/t))}=E[F(X t ,X t )]+0(n- 1 *) (8) 

Proof. For s > and y > x define h(s,x,y) = E [F(X S ,X S )\X = x,X Q = y\ . The fact 
that F(x, y) belongs to the domain of the infinitesimal generator guarantees that function 
h(s,x,y) satisfies Kolmogorov equation d s h = Ch, in particular function h(s) = h(s, 0, 0) 
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is differentiable in the variable s. Using (4) and the fact that g(n,n/t) is independent of 
X t we find 

E [F(V(n, n/t), J(n, n/t))} = E [F(X sM) , X s{n ,n/t))] 

= J E [F(X S ,X S )] P(g(n, n/t) G ds) 

R+ 

= — 7~~~~T\7 / ^(^^e-^ds. 
t»(n-l)!y v; 

R+ 



The remainder of the proof is a classical application of the stationary point method (see 
[9]). The right side of the above equation is equal to 




where we have changed the variable of integration s = t(l+u/y/n) and have used Stirling's 
formula n! = y/2nn n+ ^e~ n (l + 0(l/n)). Next, using power series expansion ln(l + x) = 
x — x 2 /2 + x 3 /3 + ... one can check that 

and combining this with (9) and the fact that h(t+tu/ 'y/n) = h(t) + h'(t)tu/\/n+o(l/ \/n) 
we obtain (8). □ 

3 Implementation 

The algorithm described in the previous section only has practical value if one is able to 
sample from the distributions of X Bl /\ and —X ei , x . It would seem that this, in itself, is 
not that much different from the problem that it purports to solve. However, it turns out 
that there are many tractable examples and in all cases this is due to the tractability of 
their Wiener-Hopf factorisations. 

Whilst several concrete cases can be handled from the class of spectrally one-sided 
Levy processes thanks to recent development in the theory of scale functions which can 
be used to described the laws of X ei /\ and — 2L ei /x ( c f- [1 1? 18]), we give here two large 
families of two sided jumping Levy processes that have pertinence to mathematical finance 
to show how the algorithm may be implemented. 
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3.1 /3-class of Levy processes 

The /3-class of Levy processes, introduced in [14], is a 10-parameter Levy process which 
has characteristic exponent 

+ ^ (b(q; 2) 1 - A 2 ) - B f a 2 + 1 - A 2 

with parameter range a, a el, Ci, c 2 , ai, a 2 , /3 2 > and Ai,A 2 G (0,3) \ {1,2}. Here 
B(x,y) = T(x)T(y)/T(x + y) is the Beta function (see [12]). The density of the Levy 
measure is given by 

Although \1/ takes a seemingly complicated form, this particular family of Levy processes 
has a number of very beneficial virtues from the point of view of mathematical finance 
which are discussed in [14]. Moreover, the large number of parameters also allows one 
to choose Levy processes within the /3-class that have paths that are both of unbounded 
variation (when at least one of the conditions a ^ 0, Ai G (2,3) or A 2 G (2,3) holds) and 
bounded variation (when all of the conditions a = 0, Ai G (0, 2) and A 2 G (0, 2) hold) as 
well as having infinite and finite activity in the jumps component (accordingly as both 
Ai, A 2 G (1, 3) or not). 

What is special about the /3-class is that all the roots of the equation A + ^(9) = 
are completly identifiable which leads to semi-explicit identities for the laws of X ei /\ and 
— 2L ei /\ as the following result lifted from [14] shows. 

Theorem 2. For A > ; all the roots of the equation 

A + ^(fl) = 

are simple and occur on the imaginary axis. They can be enumerated by {i£+ : n > 0} on 
the positive imaginary axis and '■ n > 0} on the negative imaginary axis in order of 
increasing absolute magnitude where 

C + G (0,/3 2 a 2 ), Co" e (-/3i«i,0), 

C e (ft (« 2 + n-l),(3 2 (a 2 + n)) for n > 1 

C G (/?i(-«i - n), /3i(-«i - n + 1)) for n > 1. 

Moreover, for x > 0, 



P(X ei/A G dx) = — 
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where 



1 -I 1 -I ^ 1 -I 



c = TT ' ft^±i^j and c fc - = ft( fc -^) TT ftCn-i-^O . 

n>l ' „> W fc 1 - £r 

A similar expression holds for ¥[—X ei n G dx) wit/j i/ie role of {(~ : n > 0} &emg played 
by {— Cn : n — 0} an d a i> A replaced by a 2 , /?2- 

Note that when is irregular for (0, oo) the distribution of X Bl n will have an atom 
at which can be computed from (10) and is equal to 1 — J2k>o c k- Alternatively, from 
Remark 6 in [14] this can equivalently be written as Yln>o(~(n) / Pi( n + a i)- A similar 
statement can be made concerning an atom at for the distribution of — 2L ei /x when is 
irregular for (— oo, 0). Conditions for irregularity are easy to check thanks to Bertoin [3]; 
see also the summary in Kyprianou and Loeffen [16] for other types of Ley processes that 
are popular in mathematical finance. 

By making a suitable truncation of the series (10) one may easily perform independent 
sampling from the distributions X ei /\ and 2L ei /\ as required for our Monte-Carlo methods. 



3.2 Philanthropy and General Hypergeometric Levy processes. 

The forthcoming discussion will assume familiarity with classical excursion theory of Levy 
processes for which the reader is referred to Chapter VI of [2] or Chapter 6 of [15]. 

According to Vigon's theory of philanthropy, a (killed) subordinator is called a phi- 
lanthropist if its Levy measure has a decreasing density on K + . Moreover, given any two 
subordinators Hi and H 2 which are philanthropists, providing that at least one of them 
is not killed, there exist a Levy process X such that Hi and H 2 have the same law as 
the ascending and descending ladder height processes of X, respectively. Suppose we de- 
note the killing rate, drift coefficient and Levy measures of Hi and H 2 by the respective 
triples (k,S, Ilfl-J and (A;, 5, Hh 2 )- Then [23] shows that the Levy measure of X satisfies 
the following identity 

POO ^ ^ 

^-x( x ) — / Hh 1 (x + du)Il H2 (u) + 5n Hl (x) + kIl Hl (x), x > 0, (11) 
Jo 

where ~kh x is the density of ITf^. By symmetry, an obvious analogue of (11) holds for the 
negative tail n x (x) := Ux(— oo, x), x < 0. 

A particular family of subordinators which will be of interest to us is the class of sub- 
ordinators which is found within the definition of Kuznetsov's /3-class of Levy processes. 
These processes have characteristics (c, a, ^,7) where 7 G (— 00, 0) U (0, 1), (3,c > and 
1 — a + 7 > 0. The Levy measure of such subordinators is of the type 

e a/3x 

c j^—^ l{x>0} dx - (12) 



8 



^From Proposition 9 in [14], the Laplace exponent of a /3-class subordinator satisfies 
$(0) = k + 69 + |{B (1 - a + 7 , - 7 ) - B (1 - a + 7 + 9/(3, -7) } (13) 

for 9 > where 5 is the drift coefficient and k is the killing rate. 

Let Hi and H 2 be two independent subordinators from the /3-class where for i = 1, 2, 
with respective drift coefficients ^ > 0, killing rates kj > and Levy measure parameters 
(q, «j, /3, 7j). Their respective Laplace exponents are denoted by $j, i = 1,2. In Vigon's 
theory of philanthropy it is required that kik 2 = 0. Under this assumption, let us denote 
by X the Levy process whose ascending and descending ladder height processes have the 
same law as Hi and H 2 , respectively. In other words, the Levy process whose characteristic 
exponent is given by $i(— i^)^^); 9 G R. It is important to note that the Gaussian 
component of the process X is given by 26162, see [23]. ^From (11), the Levy measure of 
X is such that 

e /3i«iu roc ^aifiiz 



1 (e^ u - 1)ti +1 J u _ x {efr* - l)^+i 



+ <52Cl (e^-l)^ +k2Cl A W u - 



dx. 



Assume first that 72 < 0, taking derivative in x and computing the resulting integrals 
with the help of [12] we find that for x > the density of the Levy measure is given by 

= -^B(p,-7 2 )e-^ (1 ^- Ql) 2 F 1 (l + 7l ,p;p-7 2 ;e-^) 



c 2 \ e ail3x d 
+ ci ( k 2 + — B(l + 72 - a 2 , -72) 7-7^ 7vu^T -5 2Ci- 



e «i/3x 



/3 v z ' ,z 7 (e^ - 1) 1+ ^ z A dx L(e^ - l) 1+ ^_ 

where p = 2 + ji + 72 — cti — a 2 - The validity of this formula is extended for 72 G (0, 1) 
by analytical continuation. The corresponding expression for x < can be obtained by 
symmetry considerations. 

We define a General Hypergeometric process to be the 13 parameter Levy process 
with characteristic exponent given in compact form 

(0) = di9 + ^a 2 9 2 + $i(-i0)$ 2 (i0), eR (14) 

where d, a G R. The inclusion of the two additional parameters d, o is largely with ap- 
plications in mathematical finance in view. Without these two additional parameters it 
is difficult to disentangle the Gaussian coefficient and the drift coefficients from param- 
eters appearing in the jump measure. Note that the Gaussian coefficient in (14) is now 
cr 2 /2 + 25 1 5 2 . The definition of General Hypergeometric Levy processes includes previously 
defined Hypergeometric Levy processes in Kyprianou et al. [17], Caballero et al. [6] and 
Lamperti-stable Levy processes in Caballero et al. [5]. 
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Just as with the case of the /3-family of Levy processes, because \l/ can be written as 
a linear combination of a quadratic form and beta functions, it turns out that one can 
identify all the roots of the equation iff (9) + A = which is again sufficient to describe 
the laws of X ei / X and — X / v 

Theorem 3. For A > ; all the roots of the equation 

\ + \ff(9) = 

are simple and occur on the imaginary axis. Moreover, they can be enumerated by {i£+ : 
n > 0} on the positive imaginary axis and {i£~ : n > 0} on the negative imaginary axis 
in order of increasing absolute magnitude where 

£o + G (0,/3(l + 72 -a 2 )),^e (-/3(l + 7 i-«!),0) 

C e (0(72 - Qt 2 + n), 0(1 + 72 - « 2 + n)) /or n > 1 
C e ( - 0(1 + 7i - «i + n), -0(7i - ai + n)) /or n > 1. 

Moreover, for x > 0, 

P(X. 1/A e ds) = - | > /it,- - ' I no) 



\fc>0 / 



where 



1+ a , € ° , S 1 + ^7 4 , i „ 1 + 



n ' /3(7i-ai+ra) , _ _ g^y-Oi+fc) TT /3(7i~a i+n) 

n>l 1 - 1 - n>l,n#* ' ~ ~ 



Moreover, a similar expression holds for P(— X Jei /\ G dx) u>z£/i £/ie ro/e o/ {£ ra : n > 0} 
replaced by {— : n > 0} and 011,72 replaced by 02,72- 

Proof. The proof is very similar to the proof of Theorem 10 in [14]. Formula (14) and 
reflection formula for the Beta function (see [12]) 

~, . . , sin(7r(;z + 7)) , . 

B(-z; -7 = B 1 + z + 7; -7 \\ " } 16 

sm(7rz) 

tell us that — )■ —00 as 9 — >■ 0(1 + 72 — 02), and since ^(0) = we conclude that 
A + (id) = has a solution on the interval 9 G (0,0(1 + 72 — 02))- Other intervals 
can be checked in a similar way (note that <f>i(z) are Laplace exponents of subordinators, 
therefore they are positive for z > 0). Next we assume that a, 81,62 > 0. Using formulas 
(14), (16) and an asymptotic result 
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which can be found in [12], we conclude that \I/(i#) has the following asymptotics as 
6 ->• +00: 

tt(i0) = -^(a 2 + 25 1 5 2 )6 2 + 0(6 1+ ^) 

/ x sin I w ( rvn -I- £ M 

^1+72 + (072+71)] 
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sin ^7r f a 2 - 72 + 



Using the above asymptotic expansion and the same technique as in the proof of Theorem 
5 in [14] we find that as n — > +00 there exists a constant C\ such that 



H = /3(n + 1 + 72 - a 2 ) + C in ^- 1 + 0{ 
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72— 1— e\ 



with a similar expression for £~. Thus we use Lemma 6 from [14] (and the same argument 
as in the proofs of Theorems 5 and 10 in [14]) to show that first there exist no other 
roots of meromorphic function A + ty(iz) except for and secondly that we have a 

factorisation 



\6 , 1 , 16 



x 1 1 ™ 1 1 + 

A 1 TT x /3(7i-ai+n) 1 TT T /3(7 2 -a 2 +n) 



n /3(7i-ai+n) 1 TT 

1 7 is x i 7 re 11 



The Wiener- Hopf factoris 4>q{0) are identified from the above equation with the help of 
analytical uniqueness result, Lemma 2 in [14]. Formula (15) is obtained from the infinite 
product representation for <f>q(9) using residue calculus. 

This ends the proof in the case cr, 61,62 > 0, in all other cases the proof is almost 
identical, except that one has to do more work to obtain asymptotics for the roots of 
A + ty(i9) = 0. We summarise all the possible asymptotics of the roots below 

= /3( n - a 2 + tu 2 ) + Cn Q2 +0(n Q2 - e ) as n -± 00. 

where the coefficients Q2 and C are presented in Table 1. Corresponding results for £~ 
can be obtained by symmetry considerations. □ 

Similar remarks to those made after Theorem 2 regarding the existence of atoms in 
the distribution of X ei /\ and — 2L ei /\ a ls° a Pply here. 

Remark 1. It is important to note that the hypergeometric Levy process is but one 
of many examples of Levy processes which may be constructed using Vigon's theory of 
philanthropy. With the current Monte-Carlo algorithm in mind, it should be possible to 
engineer other favourable Levy processes in this way. 
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Case 


U2 


C 


Q2 


rr 2 •> fl 


1 r J2 


2<5ic 2 

/3r(l+72)(<T 2 +2<5i<5 2 ) 


ru~ 1 

12 1 


a = <5i <5 9 > 


1 + To 
- 1 - ~ 12 


C2 

/3r(l+ 72 )52 




a 2 So > 5i = 


^ 1 12 


2cic 2 r(l-7i) 

i a3+ 71 -72r(i+7 2 )7io- 2 


'Yi 4- 'yo — 2 

(IT /2 ^ 


<7 2 ,<fl > 0,<5 2 = 


1 +72 


25i c 2 


72 - 1 


5 2 > 0,a = 5i = 


1+72 


C2 

/35 2 r(i+ 72 ) 


72-1 


$! > 0,<7 = 5 2 = 





sin(7T72) /3 2 7 2 (/i+d) 
7T <J 1 c 2 r(l— 72) 


-72 


a 2 > 0, ^ = 5 2 = 


1+72 


2cic 2 r(l-7i) 

/3 3+7l-72r(l+7 2 )7l cr 2 


71+72-2 


a = 5i = 62 = 


1 


^^^(k 2 + fB(l + 72 -a 2 ;- 72 )) 


-72 



Table 1: Coefficients for the asymptotic expansion of £, 



4 Extensions 

4.1 Building in arbitrary large jumps 

The starting point for the Wiener- Hopf Monte-Carlo algorithm is the distribution of X ei /\ 
and X_ ei / X , an d m section 3 we have presented two large families of Levy processes for 
which one can compute these distributions quite efficiently We have also argued the case 
that one might engineer other fit-for-purpose Wiener-Hopf factorisations using Vigon's 
theory of philanthropy. However, below, we present another alternative for extending the 
the application of the Wiener-Hopf Monte-Carlo technique to a much larger class of Levy 
processes than those for which sufficient knowledge of the Wiener-Hopf factorisation is 
known. Indeed the importance of Theorem 4 below is that we may now work with any 
Levy processes whose Levy measure can be written as a sum of a Levy measure from 
the /3-family or hypergeometric family plus any other measure with finite mass. This is 
a very general class as a little thought reveals that many Levy processes necessarily take 
this form. However there are obvious exclusions from this class, for example, cases of Levy 
processes with bounded jumps. 

Theorem 4. Let Y = {Y t : t > 0} be a sum of a Levy process X and a compound Poisson 
process such that for all t > 0, 

N t 

Y t = X t + Y,Zi 

where N = {N t : t > 0} is a Poisson process with intensity 7 and {£j : i > 1} is a sequence 
of i.i.d. random variables. Define iteratively for n > 1 

V(n,X) = V(n-1, A) + +£ n (l - 

J(n,X) = max(V(n,A) , J(n-l,X) ,V(n- 1,X) + S^\ 



12 



where V(0, A) = J(0, A) = 0, sequences {S x + 7 : n > 1} and {/}™ 7 : n > 1} are defined in 
Theorem 1, and {(3 n : n > 1} are an i.i.d. sequence of Bernoulli random variables such 
that F{(3 n = 1) = A/(7 + A). Taen 



Proof. Consider a Poisson process with arrival rate A + 7 such that points are indepen- 
dently marked with probability A/ (A + 7). Then recall that the Poisson Thinning Theorem 
tells us that the process of marked points is a Poisson process with arrival rate A. In par- 
ticular, the arrival time having index Ti is exponentially distributed with rate A. 

Suppose that T\ is the first time that an arrival occurs in the process N, in particular T\ 
is exponentially distributed with rate 7. Let e A be another independent and exponentially 
distributed random variable, and fix x G R and y > 0. Then making use of the Wiener- 
Hopf decomposition, 

(x + Y TlAex , max{y, x + Y TlAex }) 



If we momentarily set (x,y) = (V(0, A), J(0, A)) = (0,0) then by the Poisson Thin- 
ning Theorem it follows that (Y TlAex , F TlA e A ) is equal in distribution to (V(l, A), J(l, A)). 
Moreover, again by the Poisson Thinning Theorem, (Y ex ,Y ex ) is equal in distribution to 
(V(Ti, A), J{Ti, A)). This proves the theorem for the case n = 1. 

In the spirit of the proof of Theorem 1, the proof for n > 2 can be established by an 
inductive argument. Indeed, if the result is true for n = k — 1 then it is true for n = k by 
taking (x, y) = (V(k — 1, A), J{k — 1, A)) then appealing to the lack of memory property, 
stationary and independent increments of Y and the above analysis for the case that 
n — 1. The details are left to the reader. □ 

Remark 2. A particular example where the use of the above theorem is of pertinence 
is a linear Brownian motion plus an independent compound Poisson process. This would 
include for example the so-called Kou model from mathematical finance in which the 
jumps of the compound Poisson process have a two-sided exponential distribution. In 
the case that X is a linear Brownian motion the quantities X ei /\ and —X ei n are both 
exponentially distributed with easily computed rates. 

4.2 Approximate simulation of the law of (X t , X t , X_ t ) 

Next we consider the problem of sampling from the distribution of the three random 
variables (Xt, X±, X t ), which is also an important problem for applications, in particular 



(Y s{n ,x),Y g{njX) ) ± {V{T n ,X),J{T n ,X}) 



(17) 



3 

where T n = min{j > 1 : Pi = n } ■ 



i=l 
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with regard to the double-sided exit problem and, in particular, for pricing double barrier 
options. The following slight modification of the Wiener-Hopf Monte-Carlo technique 
allows us to obtain two estimates for this triple of random variables, which in many cases 
can be used to provide upper and lower bounds for certain functionals of (X t , X t , X_ t ). 

Theorem 5. Given two sequences {S^ : n > 1} and {/jj '■ n > 1} introduced in 
Theorem 1 we define iteratively for n > 1 

V(n,X) = V(n-l,X)+S^+lP 

J(n,\) = max(j(n-l,A),1/(n-l,A) + ,S A n) ) 

K(n, A) = mm(K(n- l,X),V(n,X)) (18) 

J(n, A) = max (j(n — 1, A), V(n, A) ^ 

K(n,X) = mm (^K(n - 1, X),V(n - 1, X) + I ( x n) ^j 

where V(0, A) = J(0,A) = If (0, A) = J(0, A) = K(0,A) = 0. T/jen /or any bounded 
function f(x,y, z) : 1R 3 i— )■ K which is increasing in z-variable we have 

E[f{V{n,X),J(n,X),K(n,X))) > E[/(X g(niA) ,X g(n , A) ,X g(n)A) ] (19) 
E[/(F(n,A),if(n,A),J(n,A))] < E[/(X g(n , A) , X g( „ iA) , X g(n , A) ] (20) 

Proof. From Theorem 1 we know that {V(n, A), J(n, A)) has the same distribution as 
(^g(n,A), ^g(n,x)), and, for each n > 1, K(n, A) = min{X g{fcjA) : k = 0, 1, • • • n} > X_^ n ,\)- 
The inequality in (19) now follows. The equality in (20) is the result of a similar argu- 
ment where now, for each n > 1, Kin, A) = X g , n X \ and J(n, A) = max{X g (fc )A ) : k = 
0, 1, ■ ■ ■ ,n} < X s{rhX) . □ 

Theorem 5 can be understood in the following sense. Both triples of random variables 
(V(n, A), J(n, A), K(n, A)) and (V(n, A), J(n, A), K(n, A)) can be considered as estimates 
for (X g („ jA ), X g ( n> \), 2L g r nt \)), where in the first case K(n, A) has a positive bias and in the 
second case J(n, A) has a negative bias. An example of this is handled in the next section. 



5 Numerical results 

In this section we present numerical results. We perform computations for a process X t 
in the /3-family with parameters 

(a, a, a 1 ,/3 1 ,X 1 , c 1 , a 2 , f3 2 , A 2 , c 2 ) = (a, a, 1, 1.5, 1.5, 1, 1, 1.5, 1.5, 1) 

where the linear drift a is chosen such that \l/(— i) = — r with r = 0.05, for no other 
reason that this is a risk neutral setting which makes the process {exp(X t — rt) : t > 0} 
a martingale. We are interested in two parameter sets. Set 1 has o = 0.4 and Set 2 has 
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a = 0. Note that both parameter sets give us proceses with jumps of infinite activity but 
of bounded variation, but due to the presence of Gaussian component the process X t has 
unbounded variation in the case of parameter Set 1. 

As the first example we compare computations of the joint density of (X\,Xx — X{) 
for the parameter Set 1. Our first method is based on the following Fourier inversion 
technique. As in the proof of Theorem 1, we use the fact that X ei /\ and X ei /\ — X ei n 
are independent, and the latter is equal in distribution to X ei n, to write 

P(A ei/A G dx)P(-A ei/A G dy) = P(X ei/A G dx, X ei/X - X ei/X G dy) 

= \ J e- xt F(X t G dx, X t — Xt G dy)dt 

Writing down the inverse Laplace transform we obtain 

F(X t G dx, X t - X t G dy) = ± J P(Z ei/A G dx)P(-X ei/A G dy)\~ l e Xt d\ (21) 

Ao+iK 

where Ao is any positive number. The values of analytical continuation of P(X ei / A G dx) 
for complex values of A can be computed efficiently using technique described in [14]. 
Our numerical results indicate that the integral in (21) can be computed very precisely, 
provided that we use a large number of discretization points in A space coupled with 
Filon-type method to compute this Fourier type integral. Thus first we compute the joint 
density of {X\,X\ — X\) using (21) and take it as a benchmark, which we use later to 
compare the Wiener-Hopf Monte-Carlo method and the classical Monte-Carlo approach. 
For both of these methods we fix the number of simulations M = 10 7 and the number of 
time steps iV G {20,50, 100}. For fair comparison we use 2N time steps for the classical 
Monte-Carlo, as Wiener-Hopf Monte-Carlo method with N time steps requires simulation 
of 2N random variables {S x , I A : j = 1,2, ..,N}. All the code was written in Fortran 
and the computations were performed on a standard laptop (Intel Core 2 Duo 2.5 GHz 
processor and 3 GB of RAM). 

Figure 1 presents the results of our computations. In Figure la we show our benchmark, 
a surface plot of the joint probability density function of (X 1 ,Xi — Xi) produced using 
Fourier method (21), which takes around 40-60 seconds to compute. Figures lb, lc and Id 
show the difference between the benchmark and the Wiener-Hopf Monte-Carlo result as 
the number of time steps N increases from 20 to 50 to 100. Computations take around 7 
seconds for N = 100, and 99% of this time is actually spent performing the Monte-Carlo 
algorithm, as the precomputations of the roots and the law of J A , S\ take less than 
one tenth of a second. Figure le shows the result produced by the classical Monte-Carlo 
method with N = 100 (which translates into 200 random walk steps according to our 
previous convention); this computation takes around 10-15 seconds since here we also 
need to compute the law of Xi/n, which is done using inverse Fourier transform of the 
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characteristic function of X t given in (1). Finally, Figure If shows the difference between 
the Monte-Carlo result and our benchmark. 

The results illustrate that in this particular example the Wiener-Hopf Monte-Carlo 
technique is superior to the classical Monte-Carlo approach. It gives a much more precise 
result, it requires less computational time, is more straightforward to programme and 
does not suffer from some the issues that plague the Monte-Carlo approach, such as the 
atom in distribution of X\ at zero, which is clearly visible in Figure le. 

Next we consider the problem of pricing up-and-out barrier call option with maturity 
equal to one, which is equivalent to computing the following expectation: 



Here s G [0, b] is the initial stock price. We fix the strike price K = 5, the barrier level 
b = 10. The numerical results for parameter Set 1 are presented in Figure 2. Figure 2a 
shows the graph of 7r uo (s) as a function of s produced with Fourier method similar to 
(21), which we again use as a benchmark. Figures 2b, 2c and 2d show the difference 
between the benchmark and results produced by Wiener-Hopf Monte-Carlo (blue solid 
line) and classical Monte-Carlo (red line with circles) for N G {20,50,100}. Again we 
see that Wiener-Hopf Monte-Carlo method gives a better accuracy, especially when the 
initial stock price level s is close to the barrier b, as in this case Monte-Carlo approach 
produces the atom in the distribution of X\ at zero which creates a large error. 

Figure 3 shows corresponding numerical results for parameter Set 2. In this case we 
have an interesting phenomenon of a discontinuity in 7r uo (s) at the boundary b. The dis- 
continuity should be there and occurs due to the fact that, for those particular parameter 
choices, there is irregularity of the upper half line. Irregularity of the upper half line is 
equivalent to there being an atom at zero in the distribution of X t for any t > (also at 
independent and exponentially distributed random times). We see from the results pre- 
sented in Figures 2 and 3 that Wiener-Hopf Monte-Carlo method correctly captures this 
phenomenon; the atom at zero is produced if and only if the upper half line is irregular, 
while the classical Monte-Carlo approach always generates an atom. Also, analyzing Fig- 
ures 3b, 3c and 3d we see that in this case classical Monte-Carlo algorithm is also doing 
a good job and it is hard to find a winner. This is not surprising, as in the case of param- 
eter Set 2 the process X t has bounded variation, thus the bias produced in monitoring 
for supremum only at discrete times is smaller than in the case of process of unbounded 
variation. 

Finally, we give an example of how one can use Theorem 5 to produce upper/lower 
bounds for the price of the double no-touch barrier call option 



vr uo ( s ) = e ~ r E (se 



- L {sexp(Xi)<6} 



(22) 



.dnt 



(s) = e" r E (se Xl 



K)+l 



(23) 



{sexp(Xi)<b ; sexp(.Xj)>6} 
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Function f(x, y, z) 



[se 



K)+l 



{s exp(y)<b ; s exp(z)<b} 



is increasing in both variables y 



and z, thus using Theorem 5 we find that 



,dnt/ 



7T uo (s) - e~ r E 



( se vM _ K)+l 



7rf*(s) = tt uo (s) - e" r E 



V(n,n) _ 



K) + l 



{s oxp( J(n,n))<b ; s exp(JsT(n,n))<6} 



{s cxp( J(n,n))<b ; s cxp(AT(n,n))<6} 



are the lower/upper bounds for it (s). Figure 4 illustrates this algorithm for parameter 
Set 1, the other parameters being fixed at K = 5, b = 3, b = 10 and the number 
of time steps N = 200 (400 for the classical Monte-Carlo). We see that Monte-Carlo 
approach gives a price which is almost always larger than the upper bound produced 
by the Wiener-Hopf Monte-Carlo algorithm. This is not surprising, as in the case of 
Monte-Carlo approach we would have positive (negative) bias in the estimate of infimum 
(supremum) , and given that the payoff of the double no-touch barrier option is increasing 
in infimum and decreasing in supremum this amplifies the bias. 
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(a) Fourier method benchmark (b) N = 20, WH-MC error 




(c) N = 50, WH-MC error (d) N = 100, WH-MC error 
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(a) Fourier method benchmark 



(b) N = 20, WH-MC and MC errors 
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(c) N = 50, WH-MC and MC errors (d) N = 100, WH-MC and MC errors 



Figure 2: Computing the price of up-and-out barrier option for parameter Set 1. In Figures 
(b-d) the graph of WH-MC error is blue solid line, the graph of MC error is red line with 
circles. 



21 



(a) Fourier method benchmark 



(b) N = 20, WH-MC and MC errors 




(c) N = 50, WH-MC and MC errors (d) N = 100, WH-MC and MC errors 

Figure 3: Computing the price of up-and-out barrier option for parameter Set 2. In Figures 
(b-d) the graph of WH-MC error is blue solid line, the graph of MC error is red line with 
circles. 
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Figure 4: Computing the price of the double no-touch barrier option for parameter Set 1. 
The blue solid lines represent the upper/lower bounds produced by WH-MC method, the 
red line with circles represents the MC result. 
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